import matplotlib.pyplot as plt
import matplotlib.patches as mp
import numpy as np
import qqmbr.odebook as ob
import warnings

Задача №1.

fontsize = 20
figure = plt.figure(figsize=(10, 10))
ob.axes4x4(xmin=-10.1, xmax=10, ymin=-10.1, ymax=10,fontsize=fontsize, labels=('x', 'y'))
ax = figure.gca()

fs = lambda x, y: np.array([13 * x + 6 * y, -15 * x - 5 * y]) / np.linalg.norm(np.array([13 * x + 6 * y, -15 * x - 5 * y]))
ob.mquiver(np.linspace(-10, 10, 40), np.linspace(-10, 10, 40),
           fs, color='Teal', headlength=5, scale=3, scale_units='x') 

start_points = [1, 3, 5, 7.5]
inits = [(x, 0) for x in start_points]
inits += [(-x, 0) for x in start_points]
ob.phaseportrait(lambda X: np.array([13 * X[0] + 6 * X[1], -15 * X[0] - 5 * X[1]]), 
                 inits, t=(-3, 3), n=100, linewidth=3, color='LIGHTCORAL')

ax.tick_params(axis=u'both', which=u'both',length=0)
plt.legend()
plt.show()

fontsize = 20
figure = plt.figure(figsize=(10, 10))
ob.axes4x4(xmin=-10.1, xmax=10, ymin=-10.1, ymax=10,fontsize=fontsize, labels=('x', 'y'))
ax = figure.gca()

fs = lambda x, y: np.array([4 * x - 3 * y, 3 * x + 4 * y]) / np.linalg.norm(np.array([4 * x - 3 * y, 3 * x + 4 * y]))
ob.mquiver(np.linspace(-10, 10, 40), np.linspace(-10, 10, 40),
           fs, color='Teal', headlength=5, scale=3, scale_units='x') 

start_points = [1, 3, 5, 7.5]
inits = [(x, 0) for x in start_points]
inits += [(-x, 0) for x in start_points]
ob.phaseportrait(lambda X: np.array([4 * X[0] - 4 * X[1], 3 * X[0] + 4 * X[1]]), 
                 inits, t=(-3, 3), n=100, linewidth=3, color='LIGHTCORAL')

ax.tick_params(axis=u'both', which=u'both',length=0)
plt.legend()
plt.show()

Задача №2.

Пункт (а).

fontsize = 20
figure = plt.figure(figsize=(10, 10))
ob.axes4x4(xmin=-10.1, xmax=10, ymin=-20.1, ymax=20, fontsize=fontsize, labels=('x', 'y'))
ax = figure.gca()

fs = lambda x, y: np.array([10 * x - 6 * y, 15 * x - 8 * y]) / np.linalg.norm(np.array([10 * x - 6 * y, 15 * x - 8 * y]))
ob.mquiver(np.linspace(-10, 10, 40), np.linspace(-20, 20, 40),
           fs, color='Teal', headlength=5, scale=3, scale_units='x') 

start_points = [0.1, 0.15, 0.21]
inits = [(x, 0.205) for x in start_points]
ob.phaseportrait(lambda X: np.array([10 * X[0] - 6 * X[1], 15 * X[0] - 8 * X[1]]), 
                 inits, t=(0, 5- 0.02), n=1000, linewidth=3, color='LIGHTCORAL')

ax.tick_params(axis=u'both', which=u'both',length=0)
plt.legend()
plt.show()

Пункт (b).

fontsize = 20
figure = plt.figure(figsize=(10, 10))
ob.axes4x4(xmin=-10.1, xmax=10, ymin=-10.1, ymax=10, fontsize=fontsize, labels=('x', 'y'))
ax = figure.gca()

fs = lambda x, y: np.array([-16 * x + 16 * y, -12 * x + 12 * y]) / np.linalg.norm(np.array([-16 * x + 16 * y, -12 * x + 12 * y]))
ob.mquiver(np.linspace(-10, 10, 40), np.linspace(-10, 10, 40),
           fs, color='Teal', headlength=5, scale=3, scale_units='x') 

start_points = [1, 2, 3, 4]
inits = [(x, 0) for x in start_points]
inits += [(-x, 0) for x in start_points]
ob.phaseportrait(lambda X: np.array([-16 * X[0] + 16 * X[1], -12 * X[0] + 12 * X[1]]), 
                 inits, t=(-40, 40), n=1000, linewidth=3, color='LIGHTCORAL')

ax.tick_params(axis=u'both', which=u'both',length=0)
plt.legend()
plt.show()

Пункт (c).

fontsize = 20
figure = plt.figure(figsize=(10, 10))
ob.axes4x4(xmin=-10.1, xmax=10, ymin=-10.1, ymax=10, fontsize=fontsize, labels=('x', 'y'))
ax = figure.gca()

fs = lambda x, y: np.array([8 * x - 6 * y, 4 * x - 2 * y]) / np.linalg.norm(np.array([8 * x - 6 * y, 4 * x - 2 * y]))
ob.mquiver(np.linspace(-10, 10, 40), np.linspace(-10, 10, 40),
           fs, color='Teal', headlength=5, scale=3, scale_units='x') 

start_points = [1, 3, 5, 7.5, 15]
inits = [(x, 0) for x in start_points]
inits += [(-x, 0) for x in start_points]
ob.phaseportrait(lambda X: np.array([8 * X[0] - 6 * X[1], 4 * X[0] - 2 * X[1]]), 
                 inits, t=(-40, 40), n=1000, linewidth=3, color='LIGHTCORAL')

ax.tick_params(axis=u'both', which=u'both',length=0)
plt.legend()
plt.show()

Пункт (d).

fontsize = 20
figure = plt.figure(figsize=(10, 10))
x_min = -10
x_max = 10
y_min = -10
y_max = 10
ob.axes4x4(xmin=x_min - 0.1, xmax=x_max, ymin=y_min - 0.1, ymax=y_max, fontsize=fontsize, labels=('x', 'y'))
ax = figure.gca()

fs = lambda x, y: np.array([-8 * x + 18 * y, -3 * x + 7 * y]) / np.linalg.norm(np.array([-8 * x + 18 * y, -3 * x + 7 * y]))
ob.mquiver(np.linspace(x_min, x_max, 40), np.linspace(y_min, y_max, 40),
           fs, color='Teal', headlength=5, scale=3, scale_units='x') 

start_points = [1, 5, 10, 20]
inits = [(x, 0) for x in start_points]
inits += [(-x, 0) for x in start_points]
ob.phaseportrait(lambda X: np.array([-8 * X[0] + 18 * X[1], -3 * X[0] + 7 * X[1]]), 
                 inits, t=(-10, 10), n=1000, linewidth=3, color='LIGHTCORAL')

ax.tick_params(axis=u'both', which=u'both',length=0)
plt.legend()
plt.show()

Пункт (e).

fontsize = 20
figure = plt.figure(figsize=(10, 10))
x_min = -10
x_max = 10
y_min = -10
y_max = 10
ob.axes4x4(xmin=x_min - 0.1, xmax=x_max, ymin=y_min - 0.1, ymax=y_max, fontsize=fontsize, labels=('x', 'y'))
ax = figure.gca()

fs = lambda x, y: np.array([-5 * x - 2 * y, -5 * y]) / np.linalg.norm(np.array([-5 * x - 2 * y, -5 * y]))
ob.mquiver(np.linspace(x_min, x_max, 40), np.linspace(y_min, y_max, 40),
           fs, color='Teal', headlength=5, scale=3, scale_units='x') 

start_points = [0.01, 0.8, 1.4, 2]
inits = [(x, -1) for x in start_points]
inits += [(x, 1) for x in start_points]
inits += [(-x, 1) for x in start_points]
inits += [(-x, -1) for x in start_points]
ob.phaseportrait(lambda X: np.array([-5 * X[0] - 2 * X[1], -5 * X[1]]), 
                 inits, t=(-2, 2), n=100, linewidth=3, color='LIGHTCORAL')

ax.tick_params(axis=u'both', which=u'both',length=0)
plt.legend()
plt.show()

Задача 3.

from matplotlib import animation

fontsize = 20
x_min = -10
x_max = 10
y_min = -10
y_max = 10

figure = plt.figure(figsize=(10, 10))

count = 2000
delta = 100

S_min = -20
S_max = 20
S = np.concatenate([np.linspace(S_min, S_max, count), np.array([45 / 6, 49 / 6]), np.linspace(45 / 6 - 1, 49 / 6 + 1, delta)])
S.sort()

start_points = [1, 5, 10, 20]
inits = [(x, 0) for x in start_points]
inits += [(-x, 0) for x in start_points]

def plot_phase_portrait(s):    
    plt.clf()
    ob.axes4x4(xmin=x_min - 0.1, xmax=x_max, ymin=y_min - 0.1, ymax=y_max, fontsize=fontsize, labels=('x', 'y'))
    fs = lambda x, y: np.array([5 * x - 6 * y, S[s] * x - 9 * y]) / np.linalg.norm(np.array([5 * x - 6 * y, S[s] * x - 9 * y]))
    ob.mquiver(np.linspace(x_min, x_max, 40), np.linspace(y_min, y_max, 40),
               fs, color='Teal', headlength=5, scale=3, scale_units='x') 
    ob.phaseportrait(lambda X: np.array([5 * X[0] - 6 * X[1], S[s] * X[0] - 9 * X[1]]), 
                            inits, t=(-2, 2), n=100, linewidth=3, color='LIGHTCORAL')
    plt.legend(["s = {:.3f}".format(S[s])], loc="upper left")
    return plt.plot()


anime = animation.FuncAnimation(figure, func=plot_phase_portrait, frames=len(S), blit=True)
anime.save("hw-03-3-c-long.mp4", fps=25, extra_args=["-vcodec", "libx264"])

from matplotlib import animation
from matplotlib import pyplot as plt

fontsize = 20
x_min = -10
x_max = 10
y_min = -10
y_max = 10

figure = plt.figure(figsize=(10, 10))

count = 300
delta = 100

S_min = 5
S_max = 12
S = np.concatenate([np.linspace(S_min, S_max, count), np.array([45 / 6, 49 / 6]), 
                    np.linspace(45 / 6 - 1, 49 / 6 + 1, delta)])
S.sort()
S = np.concatenate([S, S[::-1]])

start_points = [1, 5, 10, 20]
inits = [(x, 0) for x in start_points]
inits += [(-x, 0) for x in start_points]

def plot_phase_portrait(s):    
    clf()
    ob.axes4x4(xmin=x_min - 0.1, xmax=x_max, ymin=y_min - 0.1, ymax=y_max, fontsize=fontsize, labels=('x', 'y'))
    fs = lambda x, y: np.array([5 * x - 6 * y, S[s] * x - 9 * y]) / np.linalg.norm(np.array([5 * x - 6 * y, S[s] * x - 9 * y]))
    ob.mquiver(np.linspace(x_min, x_max, 40), np.linspace(y_min, y_max, 40),
               fs, color='Teal', headlength=5, scale=3, scale_units='x') 
    ob.phaseportrait(lambda X: np.array([5 * X[0] - 6 * X[1], S[s] * X[0] - 9 * X[1]]), 
                            inits, t=(-2, 2), n=100, linewidth=3, color='LIGHTCORAL')
    plt.legend(["s = {:.3f}".format(S[s])], loc="upper left")
    return plot()


anime = animation.FuncAnimation(figure, func=plot_phase_portrait, frames=len(S), blit=True)
anime.save("hw-03-3-c-short-reversed.mp4", fps=25, extra_args=["-vcodec", "libx264"])

Задача 4.

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec

w = 1
Y, X = np.mgrid[-w:w:1000j, -w:w:1000j]
U = np.sin(X) + np.exp(Y) - 1
V = np.sin(X - Y)

fontsize = 20
x_min = -1
x_max = 1
y_min = -1
y_max = 1
figure = plt.figure(figsize=(10, 10))
ob.axes4x4(xmin=x_min, xmax=x_max, ymin=y_min, ymax=y_max, fontsize=fontsize, labels=('x', 'y'))
ax = figure.gca()
ax.streamplot(X, Y, U, V, density=[1.4, 1.4], color='LIGHTCORAL', linewidth=3)

plt.tight_layout()
plt.show()

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec

w = 5
Y, X = np.mgrid[-w:w:1000j, -w:w:1000j]
U = np.sin(X) + np.exp(Y) - 1
V = np.sin(X - Y)

fontsize = 20
x_min = -5
x_max = 5
y_min = -5
y_max = 5
figure = plt.figure(figsize=(10, 10))
ob.axes4x4(xmin=x_min, xmax=x_max, ymin=y_min, ymax=y_max, fontsize=fontsize, labels=('x', 'y'))
ax = figure.gca()
ax.streamplot(X, Y, U, V, density=[1.4, 1.4], color='LIGHTCORAL', linewidth=3)

plt.tight_layout()
plt.show()

w = 1
Y, X = np.mgrid[-w:w:1000j, -w:w:1000j]
U = X + Y
V = X - Y

fontsize = 20
x_min = -1
x_max = 1
y_min = -1
y_max = 1
figure = plt.figure(figsize=(10, 10))
ob.axes4x4(xmin=x_min, xmax=x_max, ymin=y_min, ymax=y_max, fontsize=fontsize, labels=('x', 'y'))
ax = figure.gca()
ax.streamplot(X, Y, U, V, density=[1.4, 1.4], color='LIGHTCORAL', linewidth=3)

plt.tight_layout()
plt.show()

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec

w = 5
Y, X = np.mgrid[-w:w:1000j, -w:w:1000j]
U = X + Y
V = X - Y

fontsize = 20
x_min = -5
x_max = 5
y_min = -5
y_max = 5
figure = plt.figure(figsize=(10, 10))
ob.axes4x4(xmin=x_min, xmax=x_max, ymin=y_min, ymax=y_max, fontsize=fontsize, labels=('x', 'y'))
ax = figure.gca()
ax.streamplot(X, Y, U, V, density=[1.4, 1.4], color='LIGHTCORAL', linewidth=3)

plt.tight_layout()
plt.show()

def getDistance(x, y):
    return np.sqrt((x[0] - y[0]) ** 2 + (x[1] - y[1]) ** 2)


def dxdt(x, y):
    return -(np.sin(x) + np.exp(y) - 1)


def dydt(x, y):
    return -(np.sin(x - y))


def traverseToZero(point):
    dt = 1 / 10000
    
    step = 0
    
    while point[0] > 0 and point[1] > 0:
        point[0] += dxdt(point[0], point[1]) * dt
        point[1] += dydt(point[0], point[1]) * dt
        
    return point


def findSeparatrixPoint():
    steps = 10000
    
    left = 1.25
    right = 1.75
    
    minDistance = 5 + 2
    result = left
    
    step = 0
    
    for mid in np.linspace(left, right, steps): 
        start = [5., mid]
        finish = traverseToZero(start)
        distance = getDistance(finish, [0, 0])
        
        if distance < minDistance:
            minDistance = distance
            result = mid 
            
        if step % 250 == 0:
            print("step {}:\nresult = {}\ndistance = {}".format(step, result, minDistance))
        step += 1
    
    return result

print(findSeparatrixPoint())

step 0:
result = 1.25
distance = 1.7849739115506131
step 250:
result = 1.2625012501250126
distance = 1.76494237168169
step 500:
result = 1.275002500250025
distance = 1.744177205278414
step 750:
result = 1.2875037503750375
distance = 1.7224760852279404
step 1000:
result = 1.30000500050005
distance = 1.699931025793725
step 1250:
result = 1.3125062506250624
distance = 1.6764359901708235
step 1500:
result = 1.325007500750075
distance = 1.6518833824449386
step 1750:
result = 1.3375087508750876
distance = 1.6262642609800044
step 2000:
result = 1.3500100010001
distance = 1.5996693441632688
step 2250:
result = 1.3625112511251125
distance = 1.5717890965866321
step 2500:
result = 1.375012501250125
distance = 1.5426131808798733
step 2750:
result = 1.3875137513751374
distance = 1.5122311284808907
step 3000:
result = 1.40001500150015
distance = 1.480332965601611
step 3250:
result = 1.4125162516251626
distance = 1.4470085388669731
step 3500:
result = 1.425017501750175
distance = 1.411950944772639
step 3750:
result = 1.4375187518751875
distance = 1.3750557199962576
step 4000:
result = 1.4500200020002
distance = 1.3363191003351815
step 4250:
result = 1.4625212521252124
distance = 1.295352107638081
step 4500:
result = 1.475022502250225
distance = 1.2519731287600762
step 4750:
result = 1.4875237523752376
distance = 1.2059187531129612
step 5000:
result = 1.5000250025002502
distance = 1.1569452481199058
step 5250:
result = 1.5125262526252625
distance = 1.1043838374426462
step 5500:
result = 1.525027502750275
distance = 1.0478987066221388
step 5750:
result = 1.5375287528752875
distance = 0.9866200978685211
step 6000:
result = 1.5500300030003
distance = 0.9195055947978066
step 6250:
result = 1.5625312531253126
distance = 0.8451407554586309
step 6500:
result = 1.5750325032503252
distance = 0.7613272592723399
step 6750:
result = 1.5875337533753375
distance = 0.6643498573665977
step 7000:
result = 1.60003500350035
distance = 0.5468359450307857
step 7250:
result = 1.6125362536253625
distance = 0.3902361770629807
step 7500:
result = 1.625037503750375
distance = 0.01902445350647104
step 7750:
result = 1.6250875087508752
distance = 0.37154332824718456
step 8000:
result = 1.6250875087508752
distance = 0.514754749903536
step 8250:
result = 1.6250875087508752
distance = 0.620890382700418
step 8500:
result = 1.6250875087508752
distance = 0.708076318069844
step 8750:
result = 1.6250875087508752
distance = 0.7833165232895903
step 9000:
result = 1.6250875087508752
distance = 0.8500583908581011
step 9250:
result = 1.6250875087508752
distance = 0.91047237857846
step 9500:
result = 1.6250875087508752
distance = 0.9659454986708265
step 9750:
result = 1.6250875087508752
distance = 1.0173208930415119
1.6250875087508752

def traverseToZero2(point):
    dt = 1 / 100000
    
    step = 0
    
    while point[0] > 0 and point[1] > 0:
        point[0] += dxdt(point[0], point[1]) * dt
        point[1] += dydt(point[0], point[1]) * dt
        
    return point

def findSeparatrixPoint2():
    steps = 1000
    
    left = 1.622
    right = 1.626
    
    minDistance = 5 + 2
    result = left
    
    step = 0
    
    for mid in np.linspace(left, right, steps): 
        start = [5., mid]
        finish = traverseToZero2(start)
        distance = getDistance(finish, [0, 0])
        
        if distance < minDistance:
            minDistance = distance
            result = mid 
            
        if step % 100 == 0:
            print("step {}:\nresult = {}\ndistance = {}".format(step, result, minDistance))
        step += 1
    
    return result

print(findSeparatrixPoint2())

step 0:
result = 1.622
distance = 0.19331913452114777
step 100:
result = 1.6224004004004005
distance = 0.18015444960659513
step 200:
result = 1.6228008008008008
distance = 0.16592844987215152
step 300:
result = 1.6232012012012011
distance = 0.1503388701623932
step 400:
result = 1.6236016016016017
distance = 0.13290715587209628
step 500:
result = 1.624002002002002
distance = 0.11277836774958204
step 600:
result = 1.6244024024024024
distance = 0.08812344012876965
step 700:
result = 1.6248028028028028
distance = 0.052953957017226445
step 800:
result = 1.625027027027027
distance = 0.04619199612535814
step 900:
result = 1.625027027027027
distance = 0.08335707623158177
1.625027027027027

fontsize = 20
figure = plt.figure(figsize=(10, 10))
x_min = -5
x_max = 5
y_min = -5
y_max = 5
ob.axes4x4(xmin=x_min - 0.1, xmax=x_max, ymin=y_min - 0.1, ymax=y_max, fontsize=fontsize, labels=('x', 'y'))
ax = figure.gca() 

start_points = [(5, 1.6250875087508752)]
inits = start_points
ob.phaseportrait(lambda X: np.array([np.sin(X[0]) + np.exp(X[1]) - 1, np.sin(X[0] - X[1])]), 
                 inits, t=(-4.5, 100), n=1000, linewidth=3, color='Teal')

ax.tick_params(axis=u'both', which=u'both',length=0)
plt.legend()
plt.show()

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec

w = 5
Y, X = np.mgrid[-w:w:1000j, -w:w:1000j]
U = np.sin(X) + np.exp(Y) - 1
V = np.sin(X - Y)

fontsize = 20
x_min = -5
x_max = 5
y_min = -5
y_max = 5
figure = plt.figure(figsize=(10, 10))
ob.axes4x4(xmin=x_min, xmax=x_max, ymin=y_min, ymax=y_max, fontsize=fontsize, labels=('x', 'y'))

start_points = [(5, 1.6250875087508752)]
inits = start_points
ob.phaseportrait(lambda X: np.array([np.sin(X[0]) + np.exp(X[1]) - 1, np.sin(X[0] - X[1])]), 
                 inits, t=(-4.5, 100), n=1000, linewidth=3, color='Teal')

ax = figure.gca()
ax.streamplot(X, Y, U, V, density=[1.4, 1.4], color='LIGHTCORAL', linewidth=3)

plt.tight_layout()
plt.show()

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec

w = 5
Y, X = np.mgrid[-w:w:1000j, -w:w:1000j]
U = np.sin(X) + np.exp(Y) - 1
V = np.sin(X - Y)

fontsize = 20
x_min = -5
x_max = 5
y_min = -5
y_max = 5
figure = plt.figure(figsize=(10, 10))
ob.axes4x4(xmin=x_min, xmax=x_max, ymin=y_min, ymax=y_max, fontsize=fontsize, labels=('x', 'y'))

inits = [(5, 1.625)]
ob.phaseportrait(lambda X: np.array([np.sin(X[0]) + np.exp(X[1]) - 1, np.sin(X[0] - X[1])]), 
                 inits, t=(-100, 100), n=1000, linewidth=3, color='Teal')

ax = figure.gca()
ax.streamplot(X, Y, U, V, density=[1.4, 1.4], color='LIGHTCORAL', linewidth=3)

plt.tight_layout()
plt.show()

def traverse(point):
    steps = 10000000
    dt = 0.00001
    
    for _ in range(steps):
        point[0] += dxdt(point[0], point[1]) * dt
        point[1] += dydt(point[0], point[1]) * dt
        
    return point

print("Координаты точки стремления: ({}, {})".format(*traverse([5, 1.624219999])))

Координаты точки стремления: (1.0647613790473969, -2.0768312745646007)

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec

w = 6
Y, X = np.mgrid[-w:w:1000j, -w:w:1000j]
U = np.sin(X) + np.exp(Y) - 1
V = np.sin(X - Y)

fontsize = 20
x_min = -6
x_max = 5
y_min = -5
y_max = 5
figure = plt.figure(figsize=(10, 10))
ob.axes4x4(xmin=x_min, xmax=x_max, ymin=y_min, ymax=y_max, fontsize=fontsize, labels=('x', 'y'))

start_points = [(5, 1.6251)]
inits = start_points
ob.phaseportrait(lambda X: np.array([np.sin(X[0]) + np.exp(X[1]) - 1, np.sin(X[0] - X[1])]), 
                 inits, t=(-50, 100), n=1000, linewidth=3, color='Teal')

ax = figure.gca()
ax.streamplot(X, Y, U, V, density=[2, 2], color='LIGHTCORAL', linewidth=3)

plt.tight_layout()
plt.show()

print("Координаты точки стремления: ({}, {})".format(*traverse([5, 1.6251])))

Координаты точки стремления: (-5.2184239281847145, -2.0768312746171254)